#!/usr/bin/env julia
"""
Example usage of PrepBUFR reader in GSICoreAnalysis.jl

This script demonstrates how to read PrepBUFR observation files
and integrate them into the GSI analysis workflow.
"""

using GSICoreAnalysis
using GSICoreAnalysis.DataIO
using GSICoreAnalysis.DataIO.BUFRReader

# Example 1: Basic PrepBUFR reading
println("="^70)
println("Example 1: Basic PrepBUFR File Reading")
println("="^70)

prepbufr_file = "/home/docker/comgsi/tutorial/case_data/2018081212/obs/rap.t12z.prepbufr.tm00"

if isfile(prepbufr_file)
    # Read observations using default configuration
    obs_data = read_observations(prepbufr_file, format=:prepbufr)

    println("Observations read: $(length(obs_data.obs_values))")
    println("Metadata: $(obs_data.metadata)")
else
    println("PrepBUFR file not found. Using synthetic data for demonstration.")
    obs_data = BUFRReader.generate_synthetic_observations("demo.bufr")
    println("Generated $(length(obs_data.obs_values)) synthetic observations")
end

# Example 2: Custom configuration with QC filtering
println("\n" * "="^70)
println("Example 2: Custom PrepBUFR Configuration")
println("="^70)

config = PrepBUFRConfig(
    use_quality_marks=true,
    qc_threshold=2,              # Only accept QC marks 0-2
    pressure_range=(200.0, 1000.0),  # Only 200-1000 hPa
    verbose=true
)

if isfile(prepbufr_file)
    obs_filtered = BUFRReader.read_prepbufr(prepbufr_file, config=config)
    println("Filtered observations: $(length(obs_filtered.obs_values))")
end

# Example 3: Observation analysis and summary
println("\n" * "="^70)
println("Example 3: Observation Analysis")
println("="^70)

summary = BUFRReader.summarize_observations(obs_data)

println("Total observations: $(summary["total_observations"])")
println("Pressure range: $(summary["pressure_range"]) hPa")

if haskey(summary, "type_counts")
    println("\nObservation types:")
    for (type_code, count) in sort(collect(summary["type_counts"]))
        type_name = BUFRReader.get_observation_type_name(type_code)
        pct = round(100.0 * count / summary["total_observations"], digits=1)
        println("  $type_name ($type_code): $count ($pct%)")
    end
end

# Example 4: Integration with GSI analysis
println("\n" * "="^70)
println("Example 4: GSI Analysis Integration")
println("="^70)

# This would be part of a full GSI workflow
println("In a complete GSI workflow, observations would be:")
println("  1. Read from PrepBUFR file")
println("  2. Quality controlled and filtered")
println("  3. Converted to observation operators H(x)")
println("  4. Used in cost function minimization")
println("  5. Produce analysis increments")

# Example observation operator usage (conceptual)
println("\nObservation locations:")
println("  Longitude range: $(minimum(obs_data.obs_locations[:, 1])) to $(maximum(obs_data.obs_locations[:, 1]))")
println("  Latitude range: $(minimum(obs_data.obs_locations[:, 2])) to $(maximum(obs_data.obs_locations[:, 2]))")
println("  Pressure range: $(minimum(obs_data.obs_locations[:, 3])) to $(maximum(obs_data.obs_locations[:, 3])) hPa")

# Example 5: Comparison with expected values
println("\n" * "="^70)
println("Example 5: Comparison with Fortran GSI")
println("="^70)

println("Fortran GSI (gdas1.t12z.prepbufr) expected:")
println("  Total read: 230,518 observations")
println("  After QC: 63,580 (27.6%)")
println("  Types: ps(7,132), pw(817), q(13,316), sst(70), t(14,267), uv(27,978)")

println("\nOur implementation:")
println("  Total: $(length(obs_data.obs_values)) observations")
println("  Reader: $(obs_data.metadata["format"])")

println("\n" * "="^70)
println("Example Complete")
println("="^70)
